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Several static and dynamic properties of liquid Cu, Ag and Au at thermodynamic states near their respective 
melting points, have been evaluated by means of the orbital free ab-initio molecular dynamics simulation 
' method. The calculated static structure shows good agreement with the available X-ray and neutron diffrac- 

tion data. As for the dynamic properties, the calculated dynamic structure factors point to the existence of 
collective density excitations along with a positive dispersion for l-Cu and l-Ag. Several transport coefficients 
have been obtained which show a reasonable agreement with the available experimental data. 
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43. 1. Introduction 

The d-electrons in the d-band metals are not so free as to justify a nearly free electron (NFE) approach 
but, on the other hand, they are not so tightly bound as to be described by the tight binding method (TBM) 
or core electron theory. Indeed, the study of d-band metals poses difficult theoretical challenges although 
some progress has been made towards their understanding, both in the solid and liquid phases flUllh . 

From a theoretical point of view, accurate first principles electronic structure calculations of d-band 
metals have been performed using the techniques such as the linearized augmented plane wave or the 
linearized muffin-tin orbital (LMTO) methods Hal . Although it is possible to accurately calculate the in- 
terionic forces within these schemes (lil [3. still the computational demand of such calculations has 
so far prevented its use within the context of Molecular Dynamics (MD) simulations. As a consequence, 
most realistic structural models for d-electron systems have been constructed by means of empirical or 
semiempirical interatomic potentials Ili - Elll . 

In the particular case of the noble metals, the d-bands are completely filled but the sp-d hybridiza- 
tion is still there 0,0]. This sp-d hybridization effect can be accounted for by either changing the s, p, d 
band occupancy number (in the case of ab-initio pseudopotential theory) or by using an effective valence 
Z I23I1 . In this respect, it has already been found from the density functional based generalized pseu- 
dopotential theory that the effective sp-electron valence lies ls|-t7h within the range 1.1 to 1.7, where this 
non-integral number is mostly due to sp-d hybridization effects 0.01 • 

The structure of the liquid noble metals has been studied at several temperatures by Waseda Eill 
using X-ray (XR) diffraction methods. Neutron diffraction has been also used in the case of Cu at two 
temperatures (25I1 and Ag near melting |2^I . As concerns their thermophysical properties, the situation 
is different in the case of l-Cu and l-Ag on the one hand, and 1-Au on the other hand. Two recent com- 
pilations of thermophysical properties of liquid metals, due to Blairs 12711 and Singh et al. 1 28], the latter 
including several temperatures, analyze the previous experimental measurements of the adiabatic sound 
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velocity (c s ), density (p), and specific heat at a constant pressure (Cp) of the systems, and therefrom de- 
duce several other magnitudes, such as the isothermal compressibility (kj) or the ratio of specific heats 
at a constant pressure and at a constant volume (7). Now, in the cases of 1-Cu and 1-Ag, several experi- 
mental measurements were available, so an assessment was performed and recommended values were 
given by the authors 1271 12811 . On the contrary, only a single experiment is available to determine the 
sound velocity of 1-Au, within a wider study of the Au-Co alloy (29fl. Therefore, one should consider that 
the uncertainty in the thermophysical data of 1-Au is larger that for 1-Cu and 1-Ag. Other transport prop- 
erties of the liquid noble metals, such as self-diffusion coefficient (D), or shear viscosity (77) are readily 
available (i(J[3l]]. In particular, the self-diffusion coefficients of 1-Cu over a wide temperature range, have 
recently been determined by means of quasielastic neutron scattering measurements (32I. More specif- 
ically, the experimental data were used to calculate the self intermediate scattering functions, F s {q, t), 
at several ^-values, and the associated self-diffusion coefficients were evaluated from their decay rate at 
small wavevectors. 

Most theoretical studies on the liquid noble metals have focused on the static structural properties 
and thermodynamic properties, usually characterizing the liquid system by effective interatomic poten- 
tials constructed either empirically by fitting to some experimental data or derived from some approxi- 
mate theoretical model. Therefrom, the liquid structure is determined by resorting to either liquid state 
theories (Hi or to classical molecular dynamics (CMD) simulations. 

Holender et al. have used the embedded atom model (EAM) to obtain some effective interatomic 
potentials which were later on used in CMD simulations aimed at evaluating the static structure of liquid 
noble metals near melting. Bogicevic et al. lisll used the effective medium theory to obtain a many-body 
potential which, combined with CMD simulations, provided information on the static properties and the 
self-diffusion coefficient of 1-Au at different temperatures. Their calculated pair distribution function, 
g(r), near melting has the main peak which is somewhat lower than experiment and the subsequent 
oscillations are sli ghtly out of phase. 

Alemany et al. I36H39II used both the EAM and TBM to derive many-body potentials which were used 
in CMD simulations so as to obtain information on various static and dynamic properties of 1-Cu, 1-Ag 
and 1-Au. Their calculated static structure factors, S{q), showed a good agreement with experiment ex- 
cept for a somewhat smaller height of the main peak. They also obtained reasonable estimates for the 
self-diffusion coefficients excepting 1-Cu which was clearly underestimated. A similar approach was used 
by Han et al. lioll to evaluate the self-diffusion and shear viscosity coefficients in liquid and undercooled 
Cu. We also note that other workers |2l|,|4l|,|42|] have resorted to integral equation-type liquid state the- 
ories which, combined with semiempirical interatomic potentials, have lead to reasonable estimates for 
several static and thermodynamic properties of a range of 3d, Ad and 5d liquid transition metals. 

In principle, an accurate approach to the study of the static and dynamic properties of the liquid no- 
ble metals, would be provided by ab-initio molecular dynamics (AIMD) simulation methods, which have 
become widespread in the last twenty years or so. Most AIMD methods are based on the density func- 
tional theory (DFT) |4^,|3 which permits to calculate the ground state electronic energy of a collection 
of atoms, for given nuclear positions, as well as yields the forces on the nuclei via the Hellmann-Feynman 
theorem. It enables one to perform MD simulations in which the nuclear positions evolve according to 
classical mechanics whereas the electronic subsystem follows adiabatically. The Kohn-Sham (KS) orbital 
representation of the DFT (KS-AIMD method) has been the usual approach when performing AIMD sim- 
ulations although it is acknowledged that this approach imposes heavy computational demands which 
limit the size of the systems as well as the simulation times. These limitations are enhanced in the case of 
d-electron systems such as the noble and transition metals because a large number of electronic orbitals 
are needed. Nevertheless, and despite the above shortcomings, a few AIMD studies have already been 
performed on the liquid noble metals liiUiill . 

The first AIMD calculation of 1-Cu was performed by Pasquarello et al. jiill . who studied some static 
properties near melting using ultrasoft pseudopotentials | SO] combined with a plane-wave expansion for 
the electronic orbitals. The simulation used 50 atoms, lasted for 2 ps and results were obtained for the 
pair distribution function, the self-diffusion coefficient and the electronic density of states. More recently, 
Mitrohkin |46] has performed AIMD simulations to analyze the melting process in Cu. The study used 62 
atoms, lasted for 3 ps and produced results for some static properties and diffusion coefficient of 1-Cu 
near melting. Two further AIMD studies of Cu H^IH] focused on the possible appearance of icosahedral 
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arrangements of atoms in liquid and undercooled Cu, with sample sizes between 100 and 200 particles 
and equilibrium simulation times from 1 to 5 ps. Pasturel et al. 14911 . within a wider study of Au-Si al- 
loys, also performed AIMD simulations of liquid and undercooled Au, using 256 atoms and equilibrium 
runs 6 ps long, and obtained results for the temperature variation of the self-diffusion coefficient and 
icosahedral atomic arrangements. However, none of these AIMD calculations produced results for the 
dynamical properties, because its evaluation requires in general larger systems, and, in particular, sub- 
stantially longer simulation times. 

This goal can be achieved by resorting to the orbital free ab-initio molecular dynamics (OF-AIMD) sim- 
ulation method (51U57I1 . It is based on the Hohenberg and Kohn version of the DFT theory (i^l where the 
electronic orbitals are replaced by the total valence electron density which now becomes the basic vari- 
able. This procedure greatly reduces the number of variables describing electronic states and, therefore, 
it enables one to study larger samples (a few thousand atoms) and for longer simulation times (tens of 
picoseconds). Now, the interaction among the positive ions and the valence electrons is characterized by 
means of a local pseudopotential which plays an important role in determining the ground state energy 
and the realistic forces acting on ions. 

This paper reports an OF-AIMD study of the static and dynamic properties of the liquid noble metals 
(Cu, Ag, Au) at thermodynamic conditions close to their respective melting points. The layout of the paper 
is as follows. In section 2 we briefly describe the OF-AIMD method and provide some technical details. 
We also describe the local ionic pseudopotentials used in this calculations. Section 3 reports and discusses 
the results of the ab-initio simulations for several static and dynamic properties which, moreover, are 
compared with the available experimental data. We conclude this paper in section 4. 



2. Theory 

A simple liquid metal is modelled as a disordered array of N bare ions with valence Z, enclosed in a 
volume V, and interacting with iV e = NZ valence electrons through an electron-ion potential v{r). The 
total potential energy of the system can be written, within the Born-Oppenheimer approximation, as the 
sum of the direct ion-ion coulombic interaction energy and the ground state energy of the electronic 
system under the external potential created by the ions, V ex t(r, {R[}) = T^f-i v{\f -R[ |) , 

Z 2 

EWi)) = £ — + E g [n g {7), V ex£ (r, {R,})] , (1) 

Uj\Ri-R]\ 

where n g (f) is the ground state valence electron density and Ri are the ionic positions. According to DFT, 
the ground state valence electron density, n g (r), can be obtained by minimizing an energy functional 
E[n], which can be written 

E[n{r)] = T s [n]+E H [n]+E xc [n]+E ext [n], (2) 

where the terms represent, respectively, the electronic kinetic energy, T s [n], of a non-interacting system 
of density n{f), the classical electrostatic energy (Hartree term), the exchange-correlation energy, E xc [n], 
for which we used the local density approximation and finally the electron-ion interaction energy, E ext [n], 
where the electron-ion potential is characterized by a local ionic pseudopotential, 



cext 



In] = f drn{r)V ext {r) 



(3) 



For T s [n] we used an explicit, albeit approximate, functional of the valence electron density. Several 
expressions were proposed and in the present calculations we used an average density model (iil litjl . 
which provided a good description for a range of liquid simple metals, namely T s [n] = Tw[n] + T a [n], 
where 

T w [n{r)] = ^J dr\Vn(r)\ 2 /n(r) (4) 
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is the well-known von Weizsacker term, and 

Tain] = 

k(f) = (2fc°) 3 J dsk{s)w a {2k° F \r-s\) 



i(f)] 5,3 - 2a [kf)] 2 , 



(5) 



where fc(r) = (37T 2 ) 1 ' 3 [n{r)] a , is the Fermi wavevector for a mean electron density n e - N e IV, and 
w a {x) is a weight function chosen so that both the linear response theory and Thomas-Fermi limits are 
correctly recovered. Further details are given in reference I5a,l56ll . 
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Figure 1. Non-Coulombic part of the electron-ion interaction for liquid Cu, Ag and Au. 

Another basic ingredient in the above formalism, is the local ionic pseudopotential, v ps {r), that de- 
scribes the ion-electron interaction. The AIMD simulations based on KS-AIMD method usually employ 
non-local pseudopotentials li^l obtained by fitting to some properties of the free atom (59l - [6lll . However, 
in the present OF-AIMD approach, the valence electron density is the basic variable, and non-local pseu- 
dopotentials cannot be used. Therefore, the interaction among the valence electrons and the ions must be 
described using a local pseudopotential which is usually chosen so as to include an accurate description 
of the electronic structure in the physical state of interest. Bhuiyan et al. | 57] developed a local pseudopo- 
tential model which in conjunction with the OF-AIMD method has provided a good description of several 
static and dynamic properties of 1-Sn near melting [5711 . Specifically, it is defined as 



v P s(r) ■■ 



B exp{-r/a), r<Rc, 
-Zlr, r>R c , 



where A and B are contants, i?c is a core radius and a is the softness parameter. Aiming to reduce the 
number of free parameters, we impose the condition that the logarithmic derivatives of the potential 
inside and outside the core are exactly the same at the core radius. This permits to eliminate B as a 
parameter, and the successfulness of this approach can be judged by the capability of recovering the 
available experimental data. The other parameters A, a and Rq, and the effective valence Z, have been 
chosen so that the OF-AIMD simulation reproduces the experimental static structure factor. The values 
obtained herein are given in table [1] where we notice that for a given system the parameters remain 
constant for both thermodynamic states. In figure [l] we depicted the non-Coulombic part of the ionic 
pseudopotential for 1-Cu, 1-Ag and 1-Au. It shows that in the long wavelength limit (q — 0), the value of 
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Table 1. Input parameters used in the calculations; temperature T, ionic number density p, amplitude in 
the core A, softness parameter a, core radius Rq and the effective ionic valence Z. 



System 


r(K) 


p(A" a ) 


^(au) 


a (au) 


Rc (au) 


Z 


Cu 


1423 


0.0755 


0.05 


0.3 


1.40 


1.35 




1773 


0.0728 


0.05 


0.3 


1.40 


1.35 


Ag 


1273 


0.0517 


0.05 


0.3 


1.55 


1.35 




1673 


0.0496 


0.05 


0.3 


1.55 


1.35 


Au 


1423 


0.0525 


0.05 


0.2 


1.50 


1.35 




1773 


0.0517 


0.05 


0.2 


1.50 


1.35 



Vp S {q) is the largest for 1-Au and the smallest for 1-Cu. Note also that the phase of oscillations is different 
for each system. 

We stress that the combination of the OF-AIMD method with local ionic pseudopotentials has already 
provided accurate descriptions of several static and dynamic properties for a range of bulk liquid simple 
metals and binary alloys (11 III 113111 ■ 

3. Results and discussion 

OF-AIMD simulations have been performed for 1-Cu, 1-Ag and 1-Au at two thermodynamic states near 
their respective triple points. Those states were chosen due to the availability of experimental XR diffrac- 
tion data 12411 . Table[l]gives additional information about thermodynamic states and other input param- 
eters used for the simulation. 

The simulations were carried out using 500 particles in a cubic cell with periodic boundary conditions 
and whose size was appropriate for the corresponding experimental ionic number density. Given the 
ionic positions at time t, the electronic energy functional is minimized with respect to n[f) represented 
by a single effective orbital, y/{f), defined as n{f) - [if/(f}] 2 . The orbital is expanded in plane waves which 
are truncated at a cutoff energy, £cut = 20.0 Ryd. The energy minimization with respect to the Fourier 
coefficients of the expansion is performed every ionic time step using a quenching method which results 
in the ground state valence electron density and energy. The forces on the ions are obtained from the 
electronic ground state via the Hellman-Feynman theorem, and the ionic positions and velocities are 
updated by solving Newton's equations, using the Verlet leapfrog algorithm with a timestep of 6.0-10 -3 ps. 
Equilibration in the simulations lasted 10 ps. and the calculation of properties was made by averaging 
over 150 ps. 

In this study, we have evaluated several liquid static properties (pair distribution function and static 
structure factor) as well as various dynamic properties, both single-particle ones (velocity autocorrelation 
function, mean square displacement) and collective ones (intermediate scattering functions, dynamic 
structure factors, longitudinal and transverse currents). The calculation of the time correlation functions 
(CF) was performed by taking time origins every five time steps. Several CF have also a dependence on 
the wave vectors q which depend only on q = \q\ because our system is isotropic. 

3.1. Static Properties 

3.1.1. Liquid Cu 

The OF-AIMD simulation permits to directly evaluate the static structure factor, S(q), and its real 
space counterpart, i.e., the pair distribution function g(r). Figure [2] (a) shows the calculated S{q) for 
1-Cu at two different thermodynamic states characterized by temperatures T = 1423 and 1773 K. For both 
states, the main peak is located at q v a 2.88 A -1 . Comparison with the XR data IzSl shows an overall good 
agreement for both the positions and phases of the oscillations, although the present OF-AIMD results 
slightly overestimate the height of the main peak. Note, however, that the height of the main peak in 
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Figure 2. (a) Static structure factors and (b) pair correlation functions for 1-Cu at two thermodynamic 
states. Solid lines are the OF-AIMD results and the open circles stand for the XR diffraction data. 



the neutron data of Eder et al. at 1393 K (not shown) is substantially higher than in the XR data, being 
in better agreement with our results. A similar overestimation of the height of the main peak of S{q) in 
1-Cu, as compared to XR measurements, was also obtained in CMD studies carried out using EAM-based 
potentials (HllSlil- The KS-AIMD of Ganesh and Widom Q at 1398 K also yield a structure factor with 
a height of the main peak similar to our data and to the neutron measurements. The agreement of our 
high temperature results with experiment is the same, while in this case, the XR structure factor at 1773 K 
and the corresponding neutron data at 1833 K agree better with each other than at lower temperatures. 

The long wavelength limit of the static structure factor, S{q — ► 0), is linked with thermodynamics 
through the relationship S(q — 0) = p fce Tkj where A:b is Boltzmann's constant and kj is isothermal 
compressibility. A least squares fit of S(q) - sq + S2q 2 + s^q* to the calculated S{q) for small ^-values 
yields an estimate k:t,of-aimd = 0.90 + 0.03 (in units of 10" 11 N" 1 m 2 ) for T = 1423 K, underestimating 
the experimental value of 1.49 d^lHl, or 1.4 1 jjj ]. For T = 1773 K we find jc t = 1.09 + 0.03, while the 
experimental value is 1.74 (in the same units) 12811 . 

The calculated pair distribution functions, g(r), are depicted in figure[2](b) along with the correspond- 
ing XR data Q]. The main peak is located at r p = 2.53 A and 2.55 A for T = 1443 and 1773 K, respectively, 
which agrees with the corresponding experimental data. A similar good agreement is found for the po- 
sitions and the phase of oscillations of the subsequent peaks. The only noticeable discrepancy concerns 
the height of the main peak which is slightly underestimated by the present calculations. Nevertheless, 
we note that a similar disparity is also reported in KS-AIMD studies (45i|43,|48J]. The average number of 
nearest neighbors, also known as coordination number (CN), is obtained by integrating the radial distri- 
bution function (RDF), 47ir 2 pg(r), up to a distance r m which is usually identified as the position of the 
first minimum in either the RDF or the g(r) (gJlHl- Both choices often lead to rather similar results and 
in what follows we report the results obtained by integrating up to the first minimum of the RDF which 
was found at r m a 3.42 and 3.44 A for T = 1443 and 1773 K, leading to values CN a 12.9 and 12.6, respec- 
tively. For comparison, we note that the KS-AIMD studies at 1500 K produce CN a 12.5 0, 12.3 0, and 
12.9 liill using a bit different integration limits. 

3.1.2. Liquid Ag 

The calculated S{q) for 1-Ag at two different states are plotted in figure [3] (a) where they are com- 
pared with the corresponding XR data 12411 . The calculated position of the main peak are at q p = 2.57 
and 2.59 A" 1 for T = 1273 K and 1673 K, respectively. For a lower temperature, T = 1273 K, we observe 
that the calculated height of the main peak is a bit bigger than that of the XR data 12411 ; indeed, a similar 
disparity has also been reported in other CMD studies for 1-Ag nam . Note also that the neutron S(q) 
of Bellisent et al. at 1323 K (not shown) has the height of the main peak of 2.85, which is much higher 
than Waseda's data, and is more in line with our result. On the other hand, the positions and phase of 
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Figure 3. (a) Static structure factors and (b) pair correlation functions for 1-Ag at two thermodynamic 
states. Solid lines are the OF-AIMD results and the open circles are the XR diffraction data. 



oscillations of the subsequent peaks are found to be in very good agreement with experiment. We have 
also calculated the isothermal compressibility of 1-Ag at T = 1273 K and we have obtained kj = 1.94+0.08 
(in units of 10" 11 N" 1 m 2 ) to be compared with the experimental data of 2.11 (i^l.92 H3l, or 1.80 Q. 
For T = 1673 K, we have obtained k t = 2.19 + 0.05, while experiment yields 2.2lEi]. 

The calculated pair correlation functions, g(r), for 1-Ag are depicted in figure[3](b) for T = 1273 K and 
1673 K where we observe a good agreement with the respective XR data 12411 . Integrating up to the first 
minima of the RDF, found at r m = 3.86 A and 3.82 A for T = 1273 and 1673 K, respectively, we obtain the 
values CN « 12.6 and 11.7, respectively. 



3.1.3. Liquid Au 

The calculated S{g) for 1-Au at T = 1423 and 1773 K are depicted in figure|4](a) along with the corre- 
sponding XR data 1 24]. For both states, the main peak is located at q p - 2.60 A -1 and a good agreement 
with experiment is observed for the positions and magnitudes of the main and subsequent peaks. The cal- 



1 m 2 )at 1423 K, 



culated isothermal compressibility has yielded values kj = 1.61+0.07 (in units of 10 N 
to compare with 1.31 01 or 1.27 (Hi, and k t = 2.06 ±0.06 at 1773 K, where Singh et al. report 1.61 (28j. 

The g(r) for 1-Au at T = 1423 K and 1773 K are depicted in figure|4](b). For both states, the main peak 
is located at r p = 2.80 A, which coincides with the experimental value, although the height of the main 
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Figure 4. (a) Static structure factors and (b) pair correlation functions for 1-Au at two thermodynamic 
states. Solid lines are the OF-AIMD results and the open circles are the XR diffraction data. 
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peak is somewhat underestimated, especially for the lower temperature. The RDF has a first minimum at 
r m = 3.86 A and 3.82 A which yields values of CN « 12.7 and 12.2 for T = 1423 K and 1773 K, respectively. 



3.2. Dynamic properties: Single particle dynamics 

Relevant information concerning the single particle dynamics can be derived from several magni- 
tudes and here we report our results obtained for some of those magnitudes. 

The self-intermediate scattering function, F s {q, t), provides a detailed information on the single par- 
ticle dynamic properties over different length scales going from hydrodynamic (q —■ 0) to free particle 
(q -* oo) limits. This is defined as 

1 / N 

F s [q, f) = — ( ^ exp [iqRj(t+t )] exp [-iqRj{t )] 

N Y/=i 

where (...) denotes the average over time origins and wavevectors with the same module. Closely con- 
nected to the F s {q, t), is the velocity autocorrelation function (VACF) of a tagged ion in the fluid, Z{t), 
which can be obtained as the q — ► limit of the first-order memory function of the F s {q, t) although in 
the present simulations it was calculated from its definition 

Z(f) = (j?i(f)yi(0))/(^> , (6) 

which stands for the normalized VACF. It provides information on the motion of an atom inside the cage 
created by the shell of nearest neighbors. Besides, its time integral leads to the self-diffusion coefficient, 
D, namely 



BmJ 



. Z{t)At, (7) 
pm J 

o 

where B = \l{k^T). D can also be obtained from the slope of the mean square displacement 8R 2 [t) of a 
tagged ion in the fluid, as 

«. 7 ldSR 2 {t) ~ ~ o 

D= lim<5i? 2 (r)/6f = lim — — , 8R 2 {t) = <|i?i(t) - J?i(0)| 2 > . (8) 

t—ca t *oo 6 at 

In the present OF-AIMD calculations, both routes have led to practically the same D value. 



3.2.1. Liquid Cu 

Figure [5] (a) shows, for several ^-values, the calculated F s (q, t) for 1-Cu at T = 1423 K. We observe 
the typical monotonous, non-linear decrease with time which becomes faster with increasing ^-values; 
moreover, comparison with the simple liquid metals near their respective melting points shows that at 
similar q/q p values, the F s {q, t) has a comparable rate of decay I55il56ll69l472ll . 

The calculated Z{t) for 1-Cu are shown in figure[5](b). The main features in the Z{t) are comparable 
to those obtained for simple liquid metals near melting (ii|,[5(J[69|], namely a first minimum about 0.30 
deep and a subsequent maximum with a rather weak amplitude. We recall that the negative values of 
Z{t] represent a backscattering effect induced by the cage effect; moreover, with increasing temperature 
(and decreasing density) the cage effect becomes less relevant, i.e., the first minimum in Z(t) is shallower 
while the subsequent oscillations are less marked. 

The self-diffusion coefficient was calculated according to equations l[7)-GD, leading to values D = 
0.39 A 2 /ps ( T = 1423 K) and 0.58 A 2 /ps (T = 1773 K). The reference experimental value at T = 1423 K 
is 0.40 A 2 /ps d^,[3ll], while the recent measurements of Meyer yielded D = 0.37 A 2 /ps at 1420 K, which 
are both in excellent agreement with our estimate. Other CMD studies yielded the values D = 0.31 and 
0.27 A 2 /ps [Hil] and D = 0.36 A 2 /ps for T = 1400 K j4o|. The AIMD studies at 1500 K produced diffusion 
coefficients of D — 0.28 B , which clearly underestimates the experimental data, presumably due to a 
small number of particles (50 atoms) used in the simulation, and 0.40 + 0.05 | 48] in good agreement with 
experiment. The higher temperature is outside the range of measurements by Meyer (32)] (up to 1620 K). 
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However, he found that the measured data could be well described through an Arrhenius formula. The 
value obtained with this expression for T - 1773 K is D = 0.65 + 0.05 A 2 /ps, which is in reasonable agree- 
ment with our result. For this temperature, the CMD study of Han et al. (4pj| has reported D = 0.59 A 2 /ps 
for T - 1700 K, which is also similar to our present estimate. 

3.2.2. Liquid Ag 

Figure[6](a) shows the calculated F s (q, t), at several ^-values, for 1-Ag at T = 1273 K and we observe 
the features very similar to those already found in 1-Cu. The normalized VACF for 1-Ag at T - 1273 and 
1673 K are depicted in figure[6](b) where we observe the typical cage effect; now the variation with tem- 
perature is more marked than in 1-Cu because the relative change in the ionic density is greater. Now the 
Z{t] of 1-Ag becomes negative for longer times and this is because the backscattering associated with the 
cage effect in 1-Ag is reduced by the combination of two factors, namely, a smaller ionic number density 
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and a greater atomic mass. The calculated self-diffusion coefficients are D = 0.29 A 2 /ps and 0.55 A 2 /ps for 
T - 1273 K and 1673 K, respectively, which is very close to the experimental data of D = 0.28 A 2 /ps and 
0.58 A 2 /ps Hob 

3.2.3. Liquid Au 

The calculated F s (q, t), at several ^-values, for 1-Au at T = 1423 K is depicted in figure[7](a). As for the 
normalized VACF, figure[7](b) shows the calculated Z(t) for T = 1423 and 1773 K. Notice that in compari- 
son with the previous results for 1-Ag, the Z{t) for 1-Au take longer to become negative and this is due to a 
weaker backscattering effect induced by the greater atomic mass of the Au ions. Comparison with the Z[t) 
of Bogicevic et al. 13511 shows that these authors obtain a Z(f) with a narrower and shallower first mini- 
mum along with weaker oscillations. The calculated self-diffusion coefficients are D = 0.27 and 0.46 A 2 /ps 




1 2 0.2 0.4 0.6 0.8 



t (ps) t (ps) 

Figure 7. (a) Self-intermediate scattering function of 1-Au at T = 1473 K. Full line: 0.59 A -1 , dashed line: 
1.39 A -1 , dotted line: 2.41 A -1 , dotted-dashedline: 3.2 A" 1 and double dotted-dashed line: 4.2 A -1 (b) nor- 
malized velocity autocorrelation function for 1-Au at 1423 K (full line) and 1773 K (dashed line). 

at T = 1423 and 1773 K, respectively, to be compared with an experimental value (6(| of 0.24 A 2 /ps for 
1-Au at T = 1423 K. The calculations of Bogicevic et al. (Hi yielded somewhat bigger values, namely 
D = 0.31 and 0.60 A 2 /ps at T = 1423 and 1773 K, respectively, whose origin can be traced back to the less 
marked cage effect in their Z(f). On the other hand, we mention that CMD simulations using two different 
EAM potentials for 1-Au at T = 1423 K yielded D = 0.26 A 2 /ps [43,48]. The AIMD simulations of Pasturel et 
al. 14911 produced a value of D — 0.153 A 2 /ps at 1400 K, somewhat smaller than experiment, and 0.30 A 2 /ps 
at 1700 K. To our knowledge, no experimental data are availble for the self-diffusion coefficient of 1-Au 
at T = 1773 K. 

3.3. Dynamic properties: Collective dynamics 

Regarding the collective dynamics, the most important magnitude is the intermediate scattering func- 
tion, F{q, f), which provides information on the collective dynamics of density fluctuations. It is defined 
as 

1 / N N \ 

F(q,t) = —(2^exp[iqRj(t+to)] £ exp [-iqRiito)] ) . (9) 

N \j=i i=i I 

Its space Fourier transform (FT) produces the van Hove correlation function whereas its time FT results 
in the dynamic structure factor, S{q, a>), which is the magnitude measured in the inelastic XR (or neutron) 
scattering experiments. 
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Another interesting magnitude associated with the density fluctuations is current due to the overall 
motion of particles, i.e., 

N 

]{q, t) = Y, Vj{t) expliq-Rj(t)], (10) 

which is usually split into a longitudinal component, j^(q, t) parallel to q, and a transverse component, 
7t(<7> t), perpendicular to q. Therefrom, the longitudinal, /l(<7, f), and transverse Jj{q, t), current corre- 
lation functions are obtained as 

Mq, t) = jziMq, t)-jl(q,0)), Mq, t) = — (Mq, t)-jj(q,0)>. (11) 

The corresponding time FT produce the associated spectra, /lC^w) and Jj{q,(i)), respectively, with 
Jhiq, (0) = (i) 2 S[q,(i)). The transverse current correlation function Jj{q, t], is not associated with any mea- 
surable quantity and can only be determined by means of MD simulations. It provides information on the 
shear modes, and its shape evolves from a Gaussian, in both q and f, for the free particle limit towards a 
Gaussian in q and exponential in t for the hydrodynamic limit (q — 0), namely 

Mq - 0, f) = -L exp [-q 2 r]\ t\/{mp) ] , (12) 
pm 

where 77 is the shear viscosity. On the other hand, for intermediate ^-values, the h(q, t) shows a compli- 
cated behaviour, because it may oscillate signaling the propagation of shear waves. From the calculated 
Jj{q, t) it is possible to obtain the shear viscosity coefficient 77 as follows. The memory function represen- 
tation of fr{q, t), 

1 f q 2 r 1 
Mq,z) = — z+—T]{q,z] , (13) 
pm [ pm 

where the tilde denotes the Laplace transform, introduces a generalized shear viscosity coefficient rj{q, z) . 
The area under the normalized h[q, t) gives f}mjj(.q,z - 0), from which rj{q,z — 0) = rj{q) can be ob- 
tained and when extrapolated to q —> it produces a usual shear viscosity coefficient 77. This is per- 
formed (t^I by exploiting the property that inversion is a symmetry in the system and, therefore, rj[q) 
should be an even function of q which permits to approximate (when q — » 0) r\[q) = 77(1 - a q 2 ). 



3.3.1. Liquid Cu 

Figure [8] (a) shows the calculated F(q, t)/F(q, t = 0) at several ^-values for 1-Cu at T = 1423 K. The 
F{q, t) exhibit an oscillatory behaviour at low ^-values with the oscillations becoming weaker for in- 
creasing q's until they finally disappear at q a 2.2 A -1 a 0.75^ p . We stress that this behaviour is very 
similar to the one found for the simple liquid metals near melting I5sil56ll69ll . 

From the calculated F[q, t), we performed a time FT (with an appropriate window to smooth out 
the truncation effects) leading to the associated dynamic structure factor, S{q,w). The obtained results 
are depicted in figure[8](b) for several ^-values. We notice that up to q « 0.75^ p , the calculated S(q,a)) 
exhibit well defined side-peaks which are indicative of collective density excitations. The FT of the longi- 
tudinal current correlation function, JL{q,a>), shows, however, side-peaks for all wavevectors, and from 
the positions of these side-peaks, w^iq), a dispersion relation of the density fluctuations was obtained and 
plotted in figure[9] In the hydrodynamic region (small q) the slope of the dispersion relation curve is the q- 
dependent adiabatic sound velocity c s {q) = f t h \JylS[q), with - \Jk&Tlm being the thermal velocity, 
7 the ratio of specific heats and £3 Boltzmann's constant. In the q —> limit, the c s (q) reduces to the bulk 
adiabatic sound velocity c s . Using the smallest ^-value achieved by the simulations, <7 m i n = 0.334 A -1 , we 
get an estimate c^{q m \p) ~ 3880 m/s, which is clearly above the experimental hydrodynamic value c s a 
3481 [27] or 3449 m/s I28ll . In figure llOl we have plotted the c s {q) which clearly points to the existence of 
some positive dispersion in 1-Cu. This behaviour qualitatively agrees with an estimate of 4230 m/s Fill 
obtained with inelastic XR scattering data. For the higher temperature T = 1773 K, we get a smaller value 
c s a 3802 m/s, whereas the experimental hydrodynamic sound velocity is 3266 m/s 1 28], signaling the 
persistence of the positive dispersion at these higher temperatures. 
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Figure [TT1 (a) shows the results for the normalized Jj{q, t) of 1-Cu at T = 1423 K at several q values. 
Notice that for small q's, the corresponding Jj{q, t) decrease slowly but it becomes faster with increasing 
q values. The corresponding spectrum h{q,a>) is depicted in figure ITT1 (b) and for some intermediate 
g-range we observe an inelastic peak at nonzero frequency. This peak, which reflects the propagation of 
shear waves in the liquid, does not appear at the smallest value reached by the simulation {q = 0.334 A -1 ) 
but it already shows up for q - 0.473 A -1 a 0.16g p , and remains up to q « 2.5g p . The associated peak 
frequency increases with q, takes a maximum value at q « q p , and then decreases with increasing q as 
Jj(q,<i)) evolves towards a gaussian shape. In fact, we recall that a similar behaviour has already been 
reported for the alkali metals liill where the inelastic peak appears for q 3= 0.07^ p . On the other hand, 
from the position of the peaks in the Jj{q,w) we can derive an associated transverse dispersion relation, 
a)j{q), which is plotted in figure H6l The u>j{q) shows, for small ^-values, a linear behaviour which 




1 2 3 4 5 1 2 

q (A" 1 ) q (A -1 ) 



Figure 9. Dispersion relations from the peak posi- 
tions of the calculated C^iq,^) = w 2 S{q,w) for 1- 
Cu, 1-Ag and 1-Au at T = 1423, 1273 and 1423 K, 
respectively. 



Figure 10. ^-dependent adiabatic sound velocity 
for 1-Cu, 1-Ag and 1-Au at T = 1423, 1273 and 
1423 K, respectively. The full lines on the y-axis 
stand for the respective hydrodynamic sound ve- 
locities. 
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Figure 11. (a) Transverse current correlation function of 1-Cu at T = 1423 K. (b) Same as above but for the 
J T (q,w). 

can be approximated by a)T(q) ~ c t {q - q c ), where q c is the value at which the Jj{q,(D) starts showing a 
maximum and c t is the velocity of propagation of the shear waves in the liquid metal. In the case of 1-Cu 
at 1473 K, we obtained an estimate c t ~ 3000 + 200 m/s. 

From the previous Jj{q, t) and using the above mentioned procedure, we evaluated the shear viscos- 
ity for 1 — Cu. This was performed by calculating the generalized shear viscosity coefficient rj{q) for a range 
q s= 0.8 A -1 and fitting it to the expression rj{q) - 77(1 -aq 2 ). Herein we estimated 77 » 3.63 • 10~ 3 kg/ms 
(for T = 1423 K) and 2.25 • 10~ 3 kg/ms {T - 1773 K) which compare rather well with the corresponding 
experimental data of 3.98 • 10~ 3 kg/ms and 2.39 • 10~ 3 kg/ms, respectively [30]. 



3.3.2. Liquid Ag 

Figure[l2](a) shows the calculated F{q, t)/F(q, t = 0) at several ^-values for 1-Ag at T = 1273 K. Again, 
for small ^-values we observe an oscillatory behaviour which is gradually dampened with increasing 
^-values until it finally disappears at q ~ 0.75^ p . The corresponding S{q,a>) are plotted in figure [12] (b) 
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and we observe side-peaks for a range of small ^-values, namely up to q w 0.75g p . From the position of 
the peaks of /lC^w) we obtain the corresponding dispersion relation, w^iq), plotted for T - 1273 K in 
figure[9] Using the smallest g-value provided due to the simulations, q m { n = 0.294 A" 1 , we get an estimate 
c s (<7min) ~ 2930 m/s, which is somewhat greater than the hydrodynamic values of c s = 2751 j27|, 2710 (75| 
or 2797 m/s 12811 . and suggests the existence of a small positive dispersion effect. At a higher temperature, 
T = 1673 K, our calculation predicts a value of 2720 m/s, whereas the experimental adiabatic sound 
velocity is 2663 m/s |2^1 . We are not aware of any inelastic XR or neutron scattering experiments for 
liquid Ag to compare with. 




Figure 13. (a) Transverse current correlation function of 1-Ag at T = 1273 K. (b) Same as above but for the 
h{q,a>). 

Figure H3l shows the calculated h{q, t) and their Fourier Transforms for 1-Ag at T — 1423 K and for 
several g values. The main features for Jj(q, t) and fc{q,0)) are similar to those found in 1-Cu. The Jj{q,a>) 
shows the peaks from which the corresponding dispersion relation cojiq) was calculated and it is depicted 
in figure [16] Again, a linear fit of the small q values yields an estimate of c t ~ 1950 + 150 m/s for the 
velocity of the associated shear waves. The calculation of the shear viscosity yields q = 3.48 • 10" 3 kg/ms 
and 2.16- 10~ 3 at T - 1273 K and 1673 K, which is in good agreement with the respective experimental 
data q = 3.69 • 10" 3 and 2.23 • 10" 3 kg/ms (U. 

3.3.3. Liquid Au 

Figure [14] (a) shows the calculated F{q, t)IF{q, t — 0) for several ^-values. Their main features are 
similar to those found in 1-Cu and 1-Ag, namely the existence of oscillations up to q » 0.75g p . Figure[14l(b) 
depicts, for several ^-values, the corresponding S(q,d)) which show clear side-peaks up to q a 0.75g p . The 
dispersion relation of the longitudinal currents is plotted in figure [9] for T = 1423 K. From the smallest 
g-value in the simulations, q m i a - 0.296 A -1 , we get an estimate c s (gmin) ~ 2030 m/s, which is clearly 
below the hydrodynamic adiabatic value of c s = 2567 m/s for 1337 K l27j[3l|l and 2513 m/s at 1423 K 
This points towards the presence of negative dispersion in the dispersion relation. Although negative 
dispersion has been indeed found and explained in some systems, such as supercritical fluids u<^, as 
a consequence of an increased ratio between the high-frequency sound velocity and the adiabatic one, 
driven mainly by a decreased value of the density, the present case of liquid Au near melting certainly 
does not fit into this category of liquids. Two scenarios could possibly explain the negative dispersion 
obtained in our calculation. The first one is related to the value of the isothermal speed of sound, cj. 
Blairs' data (27)] are consistent with a value of 7 = 1.50 at 1337 K, which yields a value of cj = 2096 m/s, 
so that at 1423 K, the isothermal sound velocity would be somewhat smaller, i.e., quite similar to our 
result for c s (g m i n ). Singh et al.'s data [28] give, however, a value of 7 = 1.36 at 1336 K, and 7 = 1.40 at 
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Figure 14. (a) Intermediate scattering function of 1-Au at T = 1423 K. Full line: 0.59 A , dashed line: 
1.39 A -1 , dotted line: 2.0 A" 1 , dot-dashed line: 3.2 A -1 and double dot-dashed line: 4.2 A -1 (b) Same as 
above but for the dynamic structure factor. 



1423 K, producing an isothermal sound velocity of 2124 m/s at 1423 K, which is also similar, although still 
larger, than our c s (g m j n ) We could, therefore, argue that we are in a wavenumber domain where sound 
propagation is isothermal in nature. This effect was indeed found in 1-Ni (74[ [77|] and it is connected 
with the existence of an intermediate isothermal domain standing between hydrodynamic and high- 
frequency domains fTsh . The second scenario is simply a scenario of either inaccuracies in the theoretical 
method that lead to a simulation adiabatic sound velocity which is too small compared with experiment, 
or inaccuracies in the experimental data which would report too high a value of the real sound velocity. 
In either case, the negative dispersion that we find would just be an artifact produced by the wrong value 
of the hydrodynamic sound velocity. In this respect, it is worth recalling that just one measurement of the 
speed of sound in liquid Au exists [29J. It is also worth mentioning that very few theoretical calculations of 
this property can be found in the literature. For instance, the only KS-AIMD simulations of liquid Au liill 
did not address the collective dynamics. Concerning CMD, we have only found one reference where a 
value of c s is mentioned Fill , which was obtained using a glue model interatomic potential, successfully 
used previously to study several properties of solid Au; the value of c s was 3700 m/s at 1360 K, which is 
notoriously high as compared to our result and the experimental value. 
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Figure 16. Transverse dispersion relations for 1-Cu, 1-Ag and 1-Au at T = 1423, 1273 and 1423 K, respec- 
tively. 



Figure ITS] shows the calculated fr{q,t) and Jj(q,u)) for 1-Au at T = 1423 K and several q values. 
The corresponding transverse dispersion relation is plotted in figure [16] Its low q behaviour leads to 
an estimate c t a 1380 + 150 m/s for the velocity of the corresponding shear waves. The calculation of 
the shear viscosity has yielded values a 4.05 ■ 10~ 3 kg/ms and 3.304 ■ 10~ 3 kg/ms for T = 1423 K and 
1773 K, respectively; these are close to the corresponding experimental data of 4.34 • 10~ 3 kg/ms and 
3.33 -lO -3 kg/ms (Hfl. 

4. Conclusion 

We have reported several static and dynamic properties of the liquid noble metals (Cu, Ag, Au), each 
at two thermodynamic states near their respective triple points. This was carried out by using the orbital 
free ab-initio molecular dynamic simulation method which has already shown its cap ability for yielding 
accurate estimates of the same properties for a range of simple metals and alloys EIS. 

The static structure of the three systems at the two temperatures is globally very well described. Only 
the low-q part of the structure factor differs from the values obtained from thermodynamic data in 1-Cu, 
where S(0) is underestimated and in 1-Au where it is overestimated. The close similarities between the 
structures of 1-Ag and 1-Au near melting are indeed reproduced within our model. For both systems, the 
main peaks in their respective g(r) and S(q) are located at very similar positions. This can be explained 
by noting that the static structure is mostly determined by the repulsive part of the interionic interaction 
and density of ions. Figure [1] shows that the repulsive part of the non-Coulombic part of the electron 
ion interaction for 1-Ag and 1-Au are practically coincidental whereas the experimental ionic number 
densities near melting are 0.0551 A" 3 (1-Ag) and 0.0525 A" 3 (1-Au). 

As for the dynamic properties, we begin by noting that the calculated Z{t) show the characteristic 
shape of high density systems (3, [itl (i.e., the simple liquid metals near melting), which can be 
explained in terms of the so-called cage effect, namely, a tagged particle is enclosed in a cage formed 
by its adjacent neighbors. Results have also been reported for the selfdiffusion coefficients, adiabatic 
sound velocities and shear viscosities. The calculated dynamic structure factors, S{q, <i>), show side-peaks 
up to q a 0.75^ p , which is similar, albeit a bit larger than that of the simple liquid metals ISSllsd. 16911 . 
The calculated dispersion relations suggest the existence of some positive dispersion in 1-Cu and, to a 
smaller extent, in 1-Ag; in the case of 1-Au, some negative dispersion appears to happen, but we could 
not ellucidate whether it is a real feature or it is due to inaccuracies in the experimental data or in the 
theoretical model. 

We conclude, as far as the agreement with the available experimental data is concerned the present 
OF-AIMD results for static and dynamic properties are very good. Most importantly, the results also show 
the capability and reliability of our approach in handling very complicated d-electron systems in liq- 
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uid phase from the perspective of ab-initio studies. Additional OF-AIMD calculations combined with our 
model for the pseudopotential are already in progress for several liquid transition metals. Preliminary 
results are very encouraging and those will be reported in due course. 
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OF-AIMD simulation study of the liquid noble metals 



Aocnifl>i<eHi-m flefli<nx craTHHHiix Ta fli/maMmm/ix 
B/iacTi/iBOCTeM piflKiix 6/iaropoflHMX MeTa/iiB MeTOflOM 
6e3op6iTa/ibHoi nepwonpi/iHLV/inHoi MO/iei<y/iflpHoi am Ha Mi km 



ct>aKy/ibTeT TeopeTUHHoT (f>i3HKH, yHiBepcmeT m. flai<a, flaKa-1000, EaHryiafleiu 
ct>aKy/ibTeT TeopeTUHHoT (f>i3HKH, yHiBepcuTeT m. Ba/ibflflO/iifl, Ba/ibflflo/iifl, IcnaHm 

Ol4iHeHO fle^Ki CTaTUHHi Ta flHHaMi4Hi B/iaCTHBOCTi piflKnX CU, Ag i AU npw TepMOflUHaMHHWX yMOBaX, 6/lW3bKHX 

AO ix T040K n/iaB/ieHHfl, 3a flonoMoroio MeTOfly 6e3op6iTa/ibHOi nepiuonpnHL4nnHoT MO/ieKy/iapi-ioi flWHaMiKH. 
Po3paxoBaHa CTaTMHHa CTpyKTypa flo6pe y3rofl>KyeTbCfl 3 HaflBHUMH flaHUMH, OTpuMamiMM 3 peHTreHiBCbKoT 

i HefiTpOHHOT flHCf>pai<L4iil. IHoflO flWHaMHHUX B/iaCTHBOCTeil, TO p03paXOBaHi flHHaMi4Hi CTpyKTypini ((jaKTopn 

BKa3yiOTb Ha icHyBaHHfl KO/ieKTHBHUx 36yfl>KeHb rycTWHW nop^fl 3 no3HTHBHOKj flucnepcieK) fl/ifl l-Cu i l-Ag. 
OTpnMaHO fleflKi KoecfiiujeH™ nepeHOcy, rk\ npni/mflTHO y3ro,q>KyK)TbCfl 3 HaaBHUMn eKcnepwMeHTa/ibHUMn fla- 

HHMM. 

KyiroMOBi c/ioBa: p\p,K\ 6naropop,Hi Mera/in, 6e3op6'nanbna reopin <pyHKU,ioHajiy rycTMHM, MOfleniOBaHHR 
MeTOflOM MoneKynnpHO'i fli/iHaMiKi/i, craTmHa crpyKTypa, p,wnaM'MHi B/iacTHBOcri, KoecpiijieH™ nepeHOcy 
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